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1 INTRODUCTION 


Predicting genetic interactions from Boolean models of biological net¬ 
works 


Laurence Calzone,'^’^’^’^ Emmanuel Barillot,^’^’^ ® 


and Andrei Zinovyev"’^ '^’^ 


Abstract: Genetic interaction can be defined as a deviation of the phenotypic quantitative effect of a double gene mutation 
from the effect predicted from single mutations using a simple (e.g., multiplicative or linear additive) statistical model. Exper¬ 
imentally characterized genetic interaction networks in model organisms provide important insights into relationships between 
different biological functions. We describe a computational methodology allowing to systematically and quantitatively character¬ 
ize a Boolean mathematical model of a biological network in terms of genetic interactions between all loss of function and gain of 
function mutations with respect to all model phenotypes or outputs. We use the probabilistic framework defined in MaBoSS soft¬ 
ware, based on continuous time Markov chains and stochastic simulations. In addition, we suggest several computational tools 
for studying the distribution of double mutants in the space of model phenotype probabilities. We demonstrate this methodology 
on three published models for each of which we derive the genetic interaction networks and analyze their properties. We classify 
the obtained interactions according to their class of epistasis, dependence on the chosen initial conditions and phenotype. The use 
of this methodology for validating mathematical models from experimental data and designing new experiments is discussed.^ 


1 Introduction 

Genetic interaction is defined as a phenomenon by which the 
effect of a double gene mutation cannot be predicted from the 
effect of single mutations using a simple (such as additive or 
multiplicative) statistical mode l The strength of the in¬ 

teraction can be characterized by an epistatic score, which is, 
in the case of purely deleterious mutations, negative for syn¬ 
ergistic interactions (when the phenotype of a double mutation 
is significantly stronger than the expected combined effect of 
two independent single mutations), and positive for alleviating 
interactions (when the combined effect is weaker). Examples 
of synergistic interactions are synthetic lethality and synthetic 
sickness (in the case of survival-related phenotype) or synthetic 
enhancement of a phenotype^i^. An example of strong allevi¬ 
ating interaction is the suppression of an effect of one mutation 
by a second mutation (in classical genetics, such interactions 
were historically defined as “epistatic”). Genetic interactions 
in the general case of both beneficial and deleterious mutations 
can be classified into 9 groups according to various inequality 
relations between the effects of single and double mutantsi^. 

Genetic interaction networks provide important insights into 
relations between different biological functions^. Knowledge 
of genetic interactions with respect to a disease phenotype 
can provide important hints on personalized treatment strat¬ 
egy, in particular, in cance r This knowledge is cur¬ 

rently obtained by costly high-throughput screening techniques 
based on knocking-out or knocking-down genes (using siRNA 
or shRNA) in model organisms, such as yeas t worm—, 
moused or human cells^. Experimentally, one can measure 
both synthetic and synthetic dosage interactions^. Establish¬ 
ing single genetic interactions can be a result of long and te¬ 
dious work, in the case of phenotypes that are complex and 


difficult to observe such as metastasis^. 

Computational approaches have been used in order to de¬ 
rive genetic interactions from dynamical mathematical models 
or by using machine learning approaches. One of the earli¬ 
est attempts to characterize the genetic networks of the genes 
involved in metabolism was done using fiux balance analysis 
framework applied to a genome-wide reconstruction of yeast 
metabolic network—. In this work, the quantitative epistatic 
measure was introduced to characterize the genetic interactions 
as a difference between the observed effect of a double mutant 
and the multiplicative model prediction from the effect of two 
single mutation effects. It was noted that the distribution of 
the epistatic measure is tri-modal and that the interactions be¬ 
tween functional modules have a tendency for monochromatic¬ 
ity, i.e., having the same dominant sign for between-module 
interactions. In a recent paper, a similar approach was ap¬ 
plied to characterize genetic interactions with respect to mul¬ 
tiple metabolism-related phenotypes^. 

There have been many attempts to apply machine learn¬ 
ing approach for predicting genetic interactions from a subset 
of known interactions^f^. Eor instance, in yeast, the struc¬ 
ture of physical interaction networks was combined with co¬ 
expression networks; data on protein classification was used for 
predicting genetic interactions^. In worm, identical anatomi¬ 
cal expression and microarray co-expression, phenotype prox¬ 
imity, Gene Ontology annotation and presence of interlogs 
were the parameters used for fitting the logistic regression in 
order to score genetic interactions^. Decision tree-based ap¬ 
proaches trained on the structure of protein-protein interaction 
and co-expression networks in both yeast and worm were also 
used—. Short polypeptide cluster detection was utilized to pre¬ 
dict synthetic lethal interactions between genes in yeast—. Still 
in yeast, evolutionary approaches and the notion of functional 
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asymmetry allowed prediction of negative genetic interactions 
between protein complex components^ There are very few 
examples of computational predictions of genetic interactions 
in human, one of them used gene expression analysis to predict 
synthetic lethal partners of TP53 gene^. The main problem 
of most of machine learning approaches is the absence of bona 
fide negative example (absence of interaction) set for training, 
which is usually needed for a successful application of auto¬ 
mated classihcation methods^. Nevertheless, it was shown that 
machine learning methods are able to predict genetic interac¬ 
tions signihcantly better than random choice of a gene pair. 

The knowledge about molecular mechanisms involved in a 
biological phenomenon that one wishes to study can be repre¬ 
sented as a network of interacting entities—. Depending on the 
network type, the translation into a mathematical model can be 
done using an appropriate formalism (ordinary or partial dif¬ 
ferential equations, logical, rule-based modeling, etc.). These 
mathematical models can predict the effect of a perturbation, 
intrinsic or extrinsic, and anticipate the response of a drug, for 
instance. Boolean (or, more generally, logical) modeling fo¬ 
cuses on how the influences of regulatory molecules combine 
to control the expression or activity of each molecular entity 
- or process - composing the regulatory network. In a purely 
Boolean framework, each variable of the model can only take 
two values: 0 or 1 (absent/inactive or present/active). In our 
studies, we found that Boolean formalism represents a conve¬ 
nient mean of abstraction for modeling cellular biochemistry 
dynamics and verifying that the topology of the networks repre¬ 
senting the studied phenomena hts the experimentally-observed 
effects of loss or gain of function mutations on a phenotype. So 
far, there was no attempt to systematically predict genetic in¬ 
teractions using Boolean models of biological mechanisms. 

An important remark should be made with respect to any at¬ 
tempt to predict the genetic interactions computationally. Ge¬ 
netic interactions, being functional rather than physical, can 
strongly depend on the choice of both the phenotype (or model 
read-out) and the set of initial conditions used for model simu¬ 
lations. Therefore, genetic interactions can be classihed as oc¬ 
curring with respect to single versus multiple phenotypes, and 
dependent versus independent on initial conditions. With the 
mathematical model of metabolism in yeast, it was shown that 
genetic interactions synergistic with respect to one phenotype 
can become alleviating with respect to another one^. Simi¬ 
larly, depending on the set of initial conditions (accounting for 
homeostatic, physiological, nutrient-deprived, etc. conditions), 
some phenotypes represented in the model can never be reached 
or the simulations can lead to a different output with the same 
set of inputs. For example, in a model of cell fate decision pro¬ 
cess in response to TNF (or Fas) ligand activation signal, the 
cell response showed to be either survival or cell death (non- 
apoptotic and apoptotic with a higher probability for necrotic 
phenotype though) depending on the activity of some nodes of 


the model—. In a model describing the kinetics of the restriction 
point, if the G1 cell cycle phase cyclin, Cyclin D1 (CycD in the 
model), is initially active (corresponding to presence of growth 
factors), the cell enters the cycle, otherwise, it stays stuck in G1 
arrest-i^. 

In this manuscript, we suggest a quantitative methodology to 
convert a logical model of a regulatory network into a genetic 
interaction network, dehned with respect to a chosen model 
phenotype (which can be any phenotype and not only survival- 
related as it is often the case). The methodology is based on 
using the formalism of continuous time Markov chains imple¬ 
mented in MaBoSS software^^. Using published models, we 
applied our method to derive several genetic interaction net¬ 
works for the genes that compose these models. We analyze 
genetic network properties and show that they possess many 
features of experimentally-measured genetic networks. The de¬ 
rived genetic interactions reflect the functional properties of the 
mathematical models studied, so we briefly compare these pre¬ 
dicted functional relations using available databases. 

2 Methods and data 

2.1 Models used in this study 

Three published models were selected for testing the method. 
The models correspond to signalling pathways involved in can¬ 
cer with the focus on: the MAPK pathway— describing the 
crosstalk between the three mitogen-activated protein kinases: 
ERK, p38 and JNK, and their role in apoptosis and proliferation 
balance; the cell cycle with the focus on the biochemical pro¬ 
cesses regulating the restriction point-i^i^; and cell fate deci¬ 
sion between survival and death in response to extrinsic signals 
such as death receptor activation^iii. 

For each of the model, we provide the models in both GIN- 
sim— and MaBoSS— formats. Several genetic interaction net¬ 
works (GINs) per model can be constructed corresponding to 
different initial conditions and to the chosen phenotype. They 
can be found as separate Cytoscape sessions in Supplementary 
materials (Supp_Mat_GINs). 

2.2 Computing phenotype probabilities 

For each model, we computed the probability of reaching 
model phenotypes for all possible single and double mutants 
(resulting either from gain of function - modelled as fixing 
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2.3 Quantifying epistasis in double mutants 


the corresponding node value to 1, and referred to as ’’over¬ 
expression” or ”oe”, or from loss of function - modelled as fix¬ 
ing the node value to 0 and referred to as ’’deletion” or ”ko”). 
For these computations, we used both MaBoSS software and 
a set of scripts for processing the MaBoSS configuration and 
output files, implemented into BiNoM Cytoscape plugin 
MaBoSS is a C-H- software designed for simulating continu¬ 
ous/discrete time Markov processes, defined on the state tran¬ 
sition graph representing the dynamics of a Boolean network. 
MaBoSS allows the modeller to associate different rates up and 
rates down to each variable of the model when the dynamics 
is known, enabling to account for different time scales of the 
processes described by the model. Given some initial condi¬ 
tions, MaBoSS computes time trajectories by applying Monte 
Carlo kinetics algorithm (Details and examples can be found 
at: http : //maboss . curie . fr). More precisely, probabili¬ 
ties to reach a phenotype are computed as the probability for 
the variable associated to the phenotype to have the value 1, by 
simulating random walks on the probabilistic state transition 
graph. The parameters for the stochastic simulations (number 
of runs, initial conditions, maximum time, etc.) are configured 
for each simulation. The read-out can be a variable represent¬ 
ing the phenotype, a variable representing a protein or gene, or 
a combination of them. The probabilities for the selected out¬ 
puts are reported for each mutant based on predefined initial 
conditions (which can be all random). Since a state in the state 
transition graph can combine activation of several phenotype 
variables, some phenotype probabilities appear to be “mixed” 
or coupled. It is particularly the case for cyclic attractors. For 
the cell fate model, we investigated the effect of the choice of 
the initial conditions (“random” versus “physiological”) on the 
final phenotype probability distribution. The result of the sim¬ 
ulations is stored in a simple table, containing the complete set 
of mutants characterized by probabilities of all pure and mixed 
model phenotypes (in SuppJVIatJVIodels and Supp_Mat_GINs). 

2.3 Quantifying epistasis in double mutants 

2.3.1 Definition of epistasis measures. 

The results of double mutant simulations were used to quan¬ 
tify the level of epistasis between two model gene defects A 
and B with respect to a particular phenotype (j). We define the 
normalized “fitness” of a mutation (or combination of muta¬ 
tions) X with respect to a phenotype (j) as the ratio between the 
probability of the phenotype in the mutant X and the wild-type 
models. 

/f=4- (1) 

p<i> 

To fully characterize a genetic interaction, one should be able 
to characterize its strength and type. We defined the strength of 
the interaction as a deviation of the fitness of the double mutant 


from one of the four simplest statistical models frequently used 
in this context: additive, logarithmic, multiplicative and min, 
i.e., 

e^(A,B)=/f-V/(/;,/^). (2) 

where and are phenotype (j) fitness values of single gene 
defects, is the phenotype ^ fitness of the double mutant, 
and is one of the four functions: 


yf'^^^{x,y) =x + y 

[additive) 

yf^<^^ix,y)=log2ii2^-l)i2y-l) + l) 

(log) 

yf^^^{x,y) =xy 

[multiplicative) 
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To choose the best definition of \j/{x,y), the Pearson corre¬ 
lation coefficient was computed between the fitness values ob¬ 
served in all double mutants and estimated by the null model. 
The null model with maximal linear correlation was chosen: 

W{x,y) = argmaxcorr(v/^')(/^,/^),/^^), 
i=ADD,MLT,LOG,MIN. 

Note that the best definition of y/ can vary from model to 
model, from phenotype to phenotype, and even for different 
choices of initial conditions. Our simulations show that 
performs uniformly optimal or close to optimal in most of the 
simulations, having also advantage of not producing biased dis¬ 
tributions of e (see next section). 

2.3.2 Removing bias in the distribution of epistatic mea¬ 
sure values. 

After computing the distribution of epistatic measures, it can 
be observed that the peak of the distribution is shifted towards 
non-zero epistasis. This can be considered as a bias in esti¬ 
mating the null multiplicative model for quantifying the epista¬ 
sis measure (|2|i. In our experiments, it was corrected by lin¬ 
ear fitting of the observed value y = to the null model 
X = (s®® Figure [TJj). Then the epistatic measure is 

defined as: 

where a is the slope coefficient in the best linear fit estimation 
^I|y — ax\\^ —>• min. Further in the text, we refer to 
as e unless explicitly specified. 

2.3.3 Choosing the threshold for defining the set of ge¬ 
netic interactions. 

The distribution of the epistasis measures £ is asymmetric 
in many examples. Therefore, we set a threshold separately 
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Fig. 1 Illustrating epistasis measures for cell fate decision model. A) distribution of e values for three phenotypes. B) Additive model of 
epistasis, solid line shows uncorrected additive null model and dashed line shows the corrected model; an arrow shows a particular double 
mutant BAX-I-/BCL2-I-, for which the combined effect is stronger than expected by the null model (example of single-nonmonotonic genetic 
interaction, A < WT <B < AB)\ the length of the arrow equals to £{BAX + /BCL2+) in this case. C) comparison between e values for the 
case of random initial conditions and the physiological initial condition. D) comparison between e values for two different cell death 
phenotypes. 


4 




















3 RESULTS AND DISCUSSION 


for positive and negative part of epistastis measure distribu¬ 
tion (Figure [TJ\) as a multiplier of one-tailed standard devia¬ 
tions. Those genetic interactions whose strength are above k 
one-tailed standard deviations are selected, where ^ is a real 
number parameter (typically, A: = 2 as a moderately stringent 
selection criterion). 

2.3.4 Defining the type of genetic interaction. 

Since in model simulations, one can have both deleterious 
/^ < 1, neutral « 1 and beneficial /^ > 1 mutations X with 
respect to a phenotype (j), multiple possibilities arise for rela¬ 
tions between four numbers /'^, /®, and = 1 which 
cannot be simply grouped into alleviating and aggravating, as 
in the simplest case of pure deleterious mutations. We classi¬ 
fied gene interactions using the existing approach—, according 
to 75 possible inequalities between these four numbers which 
are further grouped into 9 genetic interaction classes; “sup¬ 
pressive”, “epistatic”, “conditional”, “single-nonmonotonic”, 
“additive”, “double-nonmonotonic”, “non-interactive”, “syn¬ 
thetic”, “asynthetic”. The first 4 classes in this list can be char¬ 
acterized by a direction of the interaction, i.e., mutation A is 
epistatic to B means that the effect of A completely cancels the 
effect of B (and both effects are different from the wild-type), 
and not the opposite (A B). Note that the directed genetic 
interaction maps the causal effects in opposite direction (e.g., 
mutations in downstream effectors of a phenotype can mask 
more upstream mutations). 

To define inequalities, we introduced a threshold for distin¬ 
guishing different values of fitness /, i.e., we consider two val¬ 
ues of fitness and equal, if \f^ — f^\ < where we 
typically choose 5 = 0.2. 

For example, one of the most prevalent interactions in our 
simulations is the “epistatic” (in the sense of the classical defi¬ 
nition of the notion “epistasis”) interaction which corresponds 
to inequalities B <WT < A = AB (denoting = 1 < 

fA _ or A = AB < WT < B meaning that the effects 
of single mutants are opposite with respect to the wild-type 
(one is deleterious and another is beneficial) and the effect of 
the double mutant is equal to one of the single mutants (one 
single mutant “wins”). Another interesting example is “syn¬ 
thetic” interaction type which can correspond to the inequal¬ 
ity AB < WT = A ~ B (classical “synthetic sickness”) or to 
WT =A = B< AB (“synthetic enhancement”). 

Some interaction types are counter-intuitive such as “single- 
nonmonotonic” which can correspond to the inequality A < 
WT < B < AB, when a combination of deleterious and bene¬ 
ficial mutations lead to enhancement of the phenotype stronger 
than the beneficial mutation alone. It was shown that these in¬ 
teractions are observed in real data—, and they are also ob¬ 
served in some of our simulations (see Figure|2]i. 

2.3.5 Visualizing genetic interaction network using Cy- 
toscape. 


Legend 



epistatic: A=AB<WT<B; WT<A=AB<B; WT<B<A=AB; B<WT<A=AB; 
single-nonmonotonic: A<WT<B<AB; WT<A<AB<B; AB<B<WT<A 
conditional: WT=B<A<AB; WT=B<AB<A; A<AB<WT=B; A<WT=B<AB 


suppressive: WT=A=AB<B; B<WT=A=AB 


additive: WT<B<A<AB; B<WT<AB<A; WT<A<B<AB; A<AB<WT<B; 

B<AB<WT<A;AB<B<A<WT; AB<A=B<WT 
asynthetic: WT<A=B=AB; A=B=AB<WT 
synthetic: WT=A=B<AB; WT<A=B=AB; AB<WT=A=B 
M M M M non-interactive: B=AB<WT=A; WT=A<B=AB; WT=B<A=AB; WT=A=B=AB 
double-nonmonotonic: WT<AB<B<A; AB<WT<B<A; WT=AB<B<A 


▲ 

Y 


gain-of-function 

loss-of-function 


Fig. 3 Colour code for the genetic interaction networks. The name of 
the interaction and the colour code is in accordance within. Only the 
rules found in our analyses of the three models are indicated for each 
interaction 


The selected genetic interactions are visualized in Cy- 
toscapei^ (see example with Figure 0 and Figure |2]i. The vi¬ 
sual mapping chosen distinguishes, by colour and shape, loss 
of function and gain of function single mutants. Size of the 
nodes reflects the effect on the phenotype of a single mutant, 
and the width of the edge, the epistatic effect strength of the 
corresponding double mutant. Colouring edges denotes their 
types, using the colour schema suggested befor&i^ (see Fig- 
ure[3for definition of the visualization style). 

2.3.6 Using non-linear principal component analysis for 
mapping double mutant distribution in the space of phe¬ 
notype probabilitiesThe non-linear principal manifolds were 
constructed for the distribution of all single and double mu¬ 
tants of the model in the space of computed model phenotype 
probabilities, using elastic maps metho d and ViDaExpert 
software-i^. For computation, only the mixed phenotypes with 
a probability expectation over the whole set of double mutants 
with more than 1% were selected. This results in sets of double 
mutants in multi-dimensional space for which principal mani¬ 
folds were computed (see Figure|4]i. 


3 Results and discussion 

The three Boolean models were downloaded either from The 
Cell Collective database^ or from GINsim database^. The 
stable state analysis was done in GINsim software. The mod¬ 
els were then exported in MaBoSS for simulations. Finally, 
we used some scripts embedded into BiNoM cytoscape plu¬ 
gin to automatically compute probabilities for all single and 
double mutants (including both gain of function and loss of 
function mutants for all components of each model) and vi- 
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All initial conditions 


Physiological initial conditions 
TNF=1, Fas=0 



Fig. 2 Genetic interaction networks computed for cell fate decision model, with random and physiological initial conditions and for the three 
considered phenotypes: apoptosis, necrosis and survival 
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3 RESULTS AND DISCUSSION 


3.1 Cell fate decision model 


sualize the results of paired interactions as genetic interaction 
networks. A thorough description of each model is given in 
supplementary materials (SuppMat_description_models) along 
with the Cytoscape sessions for each model, each phenotype 
and different initial conditions for one of the models (Supp- 
Mat.GINs). 


3.1 Cell fate decision model 



Fig. 4 Application of non-linear principal manifold analysis for 
visualizing the distribution of double mutants in the space of 
phenotype probabilities. The figure shows projection of phenotype 
probabilities from multi-dimensional space onto the 2D space of 
internal coordinates of the non-linear principal manifold. Each point 
corresponds to a mutant. A big violet pentagon coresponds to the 
wild-type model, triangles to single-element mutant model and 
circles to double mutants. Gradients of increase of the model 
phenotypes probabilities are shown by curved arrows. The gray color 
in the background visualizes local density of the projections onto the 
map, allowing to perform cluster analysis visually. 


Figure|2]shows the genetic interaction networks computed with 
respect to three different phenotypes (survival, apoptosis, non- 
apoptotic cell death referred to as necrosis for short )-i. 

The general shape of the epistatic measure distribution ex¬ 
hibit tri-modality (Figure[T}\), as it was previously observed in 
another modeling framework—. The cell fate decision GINs 
and the distributions of e show that the networks computed for 
different phenotypes are less similar than the networks com¬ 
puted for the same phenotype but with different initial condi¬ 
tions (Figures 12] and [T]Z!,D, with the legend for GINs given in 
Figure|3]l. Similar conclusions were made for most of the con¬ 
structed GINs in this study. 

For the physiological initial conditions with TNF=1, some 
gene alterations (and, by extension, some pathways) appear to 
be more important than when all initial conditions are consid¬ 
ered. Indeed, some of these interactions are lost in the nu¬ 
merous genetic interactions when considering all initial con¬ 
ditions. It is particularly evident for the survival phenotype. 
Overexpressing any gene from the survival pathway, which is 
described in a linear manner in this model is enough to favour 
or even force the survival phenotype. When taking in account 
all possible inputs, other pathways can help reach the survival 
phenotype: the additive effect of both RIPl and cIAP gain 
of function would be equivalent to forcing RIPlub. Single- 
nonmonotonic interactions are found numerous in the apop- 
totic and necrotic genetic networks. Unexpectedly, the gain of 
function of BCL2, which leads to a null probability of reaching 
apoptosis, together with the gain of function of BAX increases 
the apoptotic probability of BAX gain of function alone. In 
fact, BCL2 gain of function is able to block very efficiently 
both apoptosis and necrosis. If BAX gain of function promotes 
apoptosis as observed experimentally, deleting any signal from 
the necrotic (or necroptotic) pathway seems to increase apop¬ 
tosis even more. This observation confirms the mutual exclu¬ 
sive nature of the two phenotypes. In accordance with Drees 
et al.— , this type of single-nonmonotonic interactions occur 
with a high frequency in our networks but also in experimental 
data even though they are not “recognized by common genetic 
nomenclature”. 

The distribution of all single and double mutant models 
forms a set of points in the multi-dimensional space of model 
phenotype probabilities. We found it very informative to vi¬ 
sualize this set with the projection from multi-dimensional to 
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3.3 Mammalian restriction point model 


3 RESULTS AND DISCUSSION 


two-dimensional space, using advanced methods of non-linear 
data visualization such as the projection onto the principal man¬ 
ifolds constructed by the elastic maps method (Figure |4]i. In 
these visualizations, one can see that single and double mutants 
form clusters characterized by some typical phenotype proba¬ 
bility values. The cluster around the wild-type model, collects 
those mutants whose effect can be considered as neutral. Some 
clusters represent the mutants with extreme effect of induction 
of some of the phenotypes. The probability of different pheno¬ 
types changes along the non-linear directions (gradients) of in¬ 
creasing phenotype probability. Some clusters, labelled here by 
“BCL2_oe” and “ROSJco” single mutants, correspond to some 
particular states of the model (“naive survival” for “BCL2_oe” 
and a complex state combining apoptosis and “naive survival” 
for “ROSJco”: note that the last state can be artificial due to 
some irrealistic assumptions such as non-production of ROS, 
over-abundance of ATP, or impossibility of MPT). 

3.2 MAPK model 

The MAPK pathway controls several cellular processes such as 
cell cycle activation, apoptosis, survival or differentiation. The 
model of Grieco et ah— details the crosstalk between the path¬ 
ways of the three mitogen-activated protein kinases: ERK, JNK 
and p38. In response to four stimuli (EGER, EGER3, TGEbeta, 
and DNA damage), the model produces in silico the cell re¬ 
sponse in terms of proliferation, growth arrest and apoptosis 
in diverse conditions, and simulates different sets of mutations 
often found in cancer. Even though the model is generic, its 
analysis is applied to studying bladder cancerogenesis. 

Three GINs are generated using stringent conditions (inter¬ 
actions are selected above k = 3 standard deviations) for filter¬ 
ing the edges for the three phenotypes: apoptosis, growth arrest 
and proliferation. The networks are characterized by modular 
structure, in particular, for the apoptotic phenotype (Eigure |5] 
panel 1). Interestingly, interactions within some modules or 
between modules are monochromatic with respect to the type 
of the genetic interactions. Eor example, a module connect¬ 
ing several transcription factors (JUN, API, ATE2) with phos¬ 
phatase PPP2CA negatively controlling cell growth appears in 
the GINs for both the apoptosis and growth arrest phenotypes. 
All interactions inside this module are of “synthetic” type (i.e., 
synergistic). Monochromatic structure of interactions between 
modules can be seen in Eigure |5] panel 3, where the network 
can be decomposed into several modules (e.g., PTEN/p21/AKT 
versus p70/ERK/MEKO) based on the same type of interac¬ 
tions in between them. 

Genes of the apoptotic pathway such as ATM, MAX, etc. ap¬ 
pear to be hubs in the network with the emphasis on ATM and 
conditions for the two following situations: loss of function or 
gain of function of ATM and the partners that contribute to in¬ 
creasing (or compensating for the loss of) apoptosis (Eigure 


panel 1). The combination of p53 gain of function and ERK 
gain of function seems to be a good combination to improve 
the growth arrest phenotype (Eigure|5] panel 2) whereas loss of 
function of PTEN reduces the arrest caused by gain of function 
of BCL2. In the GIN for the proliferation phenotype (Eigure|5] 
panel 3), the gain of function of either MEK1_2 or ERK seems 
to be crucial in promoting proliferation, particularly in combi¬ 
nation with gain of function of AKT or loss of function of p53 
or p38, for instance. They form a hub in the network and seem 
to be very similar (symmetric) in terms of genetic interactions 
they share with the rest of the proteins of the network. 

The MAPK model is the biggest network we study here. We 
anticipate that in even bigger regulatory network models, the 
corresponding genetic interaction networks should be modular 
and provide informative hints on pathways that are activated 
with respect to a particular phenotype. Predictions about the 
co-occurrence or the mutual exclusivity between gene alter¬ 
ations could be also derived from these networks. 


3.3 Mammalian restriction point model 

This Boolean model— was adapted from a mathematical model 
based on ordinary differential equations developed by Novak 
and Tyson—. The model was built to illustrate the behaviour 
of cells exposed to cycloheximide treatments at different times 
of the cell cycle. The model describes the dynamics of the re¬ 
striction point situated in late G1 after which the cell commits 
to division even if treated by the drug. 

Eor this small model, the GINs are easier to interpret biologi¬ 
cally (Eigure|6]l. The model is built such that if the cell does not 
receive any external growth signals, of which CycD is the sen¬ 
sor, it remains stuck in G1 cell cycle phase. Therefore, neither 
CycD nor Rb are included in these networks as their gain or loss 
of function would automatically lead to forcing or deleting the 
phenotypes. The gain of function of the cell cycle inhibitor p27 
is counteracted by the gain of function of downstream cyclins 
such as CycA and CycE. Similarly, if both inhibitors of the G2 
and M cyclins are deleted, Cdc20 and cdhl, it is equivalent to 
overexpressing the cyclins and the cells can no longer arrest. 
A similar mechanism is achieved by overexpressing E2E and 
deleting cdhl. The role of cdhl seems to be more prevalent in 
degrading the cyclins. Note that cdhl and Cdc20 are in both 
genetic interaction networks for growth arrest and proliferation 
because the two read-outs are symmetric. The loss of function 
of both Cdc20 and cdhl leads to a very low probability of ar¬ 
resting the cycle, and a very high probability for proliferating. 
The two phenotypes are mutually exclusive. 
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Fig. 5 Genetic interaction networks computed for MAPK model, 
with random initial conditions and for the three phenotypes: 
apoptosis, growth arrest and proliferation 



Fig. 6 Genetic interaction networks computed for restriction point 
model, with random initial conditions and for the two phenotypes: 
growth arrest and proliferation 
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4 CONCLUSIONS 


3.4 Comparison with experimentally derived genetic in¬ 
teractions 

We performed two types of comparisons; first, we compared 
the genetic interactions from our method to available experi¬ 
mental results, and second, we compared the genetic interac¬ 
tions between models. 

We have compared the results from each of the examples we 
have chosen in this analysis with genetic interactions listed in 
BioGRID database^. In the database, we queried the genes 
that appeared as participating in pairs of genetic interactions 
in a significant manner in our three models. We found that in 
the MARK model, TP53 and MDM2 interactions came out in 
both BioGRID and our study: TPS3 and MDM2 were identi¬ 
fied in a phenotypic suppression type of genetic interaction in 
BioGRID and we showed that overexpression of both TP53 and 
of MDM2 led to a suppressive genetic interaction with respect 
to the apoptosis phenotype. The pair ATM and TP53 seems to 
be involved in a phenotypic enhancement in BioGRID, but was 
not found in our study. In the cell fate model, we listed three 
phenotypic suppressions between XIAP and CASP8, IKKl and 
TNF, and BCL2 and CASP8. The first two were confirmed in 
our analysis: overexpression of XIAP and of CASP8 lead to 
an epistatic interaction with respect to apoptosis in the TNF- 
activated signal, and deletion of IKKl and deletion of TNF 
lead to an epistatic interaction with respect to the necrosis 
(NonACD) phenotype in the TNF-activated signal. Also, over¬ 
expression of IKKl and deletion of TNF lead to an epistatic 
interaction with respect to the survival phenotype in the TNF- 
activated signal. The last interaction was not identified with our 
method. In the mammalian restriction point model, there was 
only one interaction that appeared in BioGRID and involved 
a phenotypic enhancement between p21 and p27 which was 
not found in our analysis. More details can be found in sup¬ 
plementary materials, SuppMat_Analysis_BioGRID. In conclu¬ 
sion, the comparison showed that some interactions predicted 
by our method were indeed confirmed in BioGRID database. 
This type of comparison can serve to validate Boolean models 
developed for various molecular mechanisms with respect to 
known genetic interactions and provide additional constraints 
on the choice of model network topology, logical rules and rate 
parameters. Of course, in this analysis one should take into ac¬ 
count incompleteness of our knowledge on genetic interactions. 

We also compared more in detail the results of the genetic 
interactions among the three examples. Unfortunately, there 
was no overlap between the three models since the only com¬ 
mon gene was BCL2 between the cell fate and the MAPK 
models. We then looked more carefully at the genetic inter¬ 
actions between phenotypes but for each model individually. 
With this comparison, we identified the complementary role of 
some genes in the networks and confirmed findings from the 
initial publications. The results can be found in supplementary 


materials, SuppMat_comparison_phenotypes. 

4 Conclusions 

In this manuscript, we suggest a methodology for converting 
a logical mathematical model with a set of initial conditions 
into the corresponding genetic interaction network characteriz¬ 
ing the behaviour of all single and double mutants in terms of 
phenotype probabilities. The advantage of the methodology is 
in that it allows; 

1) estimating and classifying possible functional interactions 
between the different elements composing the model; 

2) distinguishing extreme cases of mutations amplifying or 
masking each other and, based on this, suggesting intervention 
points in order to achieve a desired phenotype (such as in—); 

3) suggesting experimental designs from the logical models; 

4) detecting controversial (non-intuitive) properties of mu¬ 
tants with respect to expected phenotypes such as nonmono¬ 
tonic genetic interactions; 

5) comparing quickly similar logical network models in 
terms of their functional properties; 

6) validating the model and comparing different models us¬ 
ing available screenings for genetic interactions (such as syn¬ 
thetic lethality screens). 

The last point deserves further development. We aim at ex¬ 
tending our methodology using existing databases containing 
genetic interactions (similar to what we did with BioGRID^) 
for matching the model predictions with genetic interactions or 
single mutation phenotypes known from the literature or from 
screenings. Moreover, similar to the methodology of param¬ 
eter fitting in constructing chemical kinetic models, one can 
fit the kinetic rates defined in our continuous-time discrete ap¬ 
proach— in order to optimize the set of model predictions. 
Another set of experimental data that could be used with this 
approach is high-throughput cancer data, such as large-scale 
mutation landscapes that are collected for series of tumours. 
Patterns of co-occurrence or mutual exclusivity of mutations 
can reflect action of genetic interactions in cancer cells. For 
example, synthetically lethal interactions can lead to the pat¬ 
tern of mutual exclusivity since cancer cells possessing both 
synthetically lethal mutations will be eliminated from the cell 
population. Using these data for interpretation and validation 
of model-based predictions requires the development of a sta¬ 
tistical methodology for detecting statistical patterns in high- 
throughput data. 

Genetic interaction networks reconstructed from log¬ 
ical mathematical models possess many properties of 
experimentally-measured networks. They are characterized by 
a variety of types of genetic interactions (with predominance 
of masking, e.g., epistatic interactions), modular structure for 
sufficiently big discrete models (Figure H), with some modules 
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characterized by monochromaticity for within-module interac¬ 
tions as well as between-module interactions. Sets of genetic 
interactions are highly dependent on the phenotype with re¬ 
spect to which they are defined and, to less extent, sensitive 
to the initial conditions (in other words, to the molecular con¬ 
text) chosen for performing simulations. These properties make 
the obtained genetic interaction networks a good model for the 
experimentally-measured ones. 

Therefore, we believe that the suggested methodology will 
contribute to the toolbox of computational approaches in sys¬ 
tems biology, connected to mathematical modeling of cellular 
mechanisms. 
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